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Abstract 

The problem of structured noise suppression is addressed by i)modelling the subspaces 
hosting the components of the signal conveying the information and ii)applying a non- 
extensive nonlinear technique for effecting the right separation. Although the approach is 
applicable to all situations satisfying the hypothesis of the proposed framework, this work 
is motivated by a particular scenario, namely, the cancellation of low frequency noise in 
broadband seismic signals. 



1 Introduction 

The problem of structured noise suppression concerns the elimination of signal components 
produced by phenomena interfering with the observations of interest. This problem can be 
addressed by linear techniques provided that the subspaces hosting the signal components are 
complementary and well separated [1-4]. More precisely, if a signal represented by the ket 
I/) is produced as the superposition of two components G Si and I/2) G ^2, provided 
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that 5i n 1S2 = {0}, the components of the superposition |/) = + I/2) can be separated 
by an obhque projection. Even when this condition is theoretically fulfilled, if the subspaces 
<Si and S2 are not well separated, the concomitant linear problem for extracting one of the 
signal components may be ill posed, which causes the failure to correctly split the signal by a 
linear operation. Hence, nonlinear techniques for determining a subspace V C iSi, such that 
1/1) e V, and the projection onto V along ^2 is well posed, have been considered [3-5]. In those 
publications the theoretically complementary subspaces »Si and S2 are assumed to be known. 
Nevertheless, the condition Si f] S2 — {0} is strong and the possibihty of meeting it depends 
on the ability to generate the right model for the subspaces. Unfortunately, the modelling of 
the complementary subspaces by pure physical considerations is not always possible and one 
needs to relay on more general mathematical modelling. 

Although the technique for subspace modelling we introduce here is applicable to different 
situations, the work is motivated by a particular problem relevant to the processing of seismic 
signals. In the nearshore these signals may be affected by a population of low-frcqucncy waves 
called infragravity waves [6]. This type of noise may be also unavoidable in bottom broadband 
seismic observations [7]. The interested reader is refereed to [8] for explanations on how in- 
fragravity waves are generated. We restrict our consideration to the problem of reducing that 
type of structured low frequency noise from broadband seismic signals. 

Our purpose is twofold. We aim at i) mathematically modelling the subspaces to represent 
the signal components ii)provide a sparse enough representation of the signals so as to make 
sure that the correct splitting can be realized. 

Under the hypothesis that one of the signals components lies in the subspace of low fre- 
quency signals, we determine the subspace of the other component in an adaptive manner. We 
assume that such a component belongs to an unknown spline space and determine the knots 
characterizing the space by taking into account the curvature points of the signal in hand. In 
that sense, the space is 'adapted' to the particular signal being analyzed. In line with [4] we 
tackle the problem of finding the representation of this component through the minimization 
of the g-norm like quantity, which is closely related to the non-extensive entropy introduced as 
ingredient of a thermodynamic framework in the seminal paper by Tsallis [9] and ever since 
broadly applied in physics [9-16] and other disciplines [17]. 

The paper is organized as follows: In Section 2 we address the problem of subspaces mod- 
elling and discuss the non-extensive nonlinear technique yielding the right signal splinting. A 
numerical simulation concerning the filtering of low frequency noise from a seismic signal is 
presented in Section 3. The conclusions are drawn in Section 4. 

2 Adaptive subspace modelling for structured noise fil- 
tering 

As already mentioned, we are concerned with the problem of separating from a signal those 
components which are not relevant to the phenomenon of interest. For simplicity we consider 
that a signal |/) is the superposition of only two components and the goal is to find a suitable 
model for the subspaces hosting such component. Since out work is motivated by the specific 
problem of filtering low frequency noise from a seismic signal, we further assume that the 
subspace representing that type of structured noise is spanned by a few Fourier functions. In 
accordance with previous works we denote such a subspace as and consider it to be fixed. 
The composed signal is the superposition |/) = |/v) + |/w-l) with |/y) e V. The goal is to 
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model the subspace V fulfilling the theoretical condition V fl W"*" = {0}, regardless of the fact 
that the two subspaces may be too close to each other for the signal separation to be obtained 
via a linear approach. We allow for this difficulty by introducing the additional hypothesis that 
l/v) is well represented in a subspace of V, which is tantamount to assuming that |/v) has a 
'sparse enough' representation in V. We further assume that |/v) is well approximated in a 
dedicated spline space to be constructed as described in the next section. 



2.1 Finding the appropriate spline space 

Let us start by stating the few definitions on spline spaces which are needed for setting up 
our mathematical framework. For a complete treatment of splines we refer to the fundamental 
books [18-20]. 

Definition 1. Given a finite closed interval [c, d] we define a partition of [c, d] as the finite set 
of points 

A := {xjj^Q^, e N, such that c — xq < xi < ■ ■ ■ < xn < xn+i — d. (1) 

We further define N subintervals Ii,i — 0,...,N as: li — [xj, Xj+i), i = 0, . . . , N — 1 and 
In 

Definition 2. Let be the space of polynomials of degree smaller or equal to m & Nq — 
N U {0}. Let m be a positive integer and define 

SU^) = {feC"^-'[c,d] : /|;^Gn„_i,z = 0,...,iV}, (2) 

where /|/. indicates the restriction of the function f on the interval li. 

An extended partition with single inner knots associated with 5'^(A) is a set A = {yi}^^^^ 
such that 

Vm+i = Xi, i=l,...,N, Xi < ■ ■ ■ < Xn 

and the first and last m points yi < ■ ■ ■ < y-m ^ c, d < ym+N+i < ■ ■ ■ < y2m+N can be 
arbitrarily chosen. With each fixed extended partition A there is associated a unique B-spline 
basis for Sm{^), that we denote as {BmjjJJi^ ■ The B-sphne J can be defined by the 
recursive formulae [18]: 



tj ^ X *\ tj-^\^ 




Bm,j{x) = ^ Bm-l,j{x) + ^ Bm-l,j+l{x). 

yj+m—l Vj yj+m Vj+l 

For each order, m, the corresponding spline space is determined by the number and position of 
the knots. In Fig 1 we show B-spline basis for two different cubic sphne spaces (m = 4) having 
the same number of knots but located at different positions. 

In order to find the appropriate partition A giving rise to the appropriate spline space ^'^(A) 
to represent a given signal /, we first determine the critical points of the signal's curvature, 
i.e., we find the set T defined as 
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The entries in T are chosen as the initial knots of A. Extra knots are obtained by subdivisions 
between consecutive knots in T so as to generate a partition A with the desired number of 
knots. An algorithm for implementing this procedure on a signal given as a discrete piece of 
data is outlined in [21]. We would like to be able to use this procedure on the signal we need 
to represent, namely the component |/v), but, of course we do not have access to this signal; 
our goal is to find it!. Thus, in line with [3,4] we proceed as explained below. 

Since in our framework W"*" is fixed and known, we can construct the orthogonal projector 
onto W = (W"'")"'" that we denote P)^. Assuming now that V = S'm(A), for some order m and 
some partition A, we can use B-splines to span the space so that for |/v) G S'.m(A) we have 

M M 
fv{x) = (x|/v) = ^Ci(x|5i_^) = ^^dBi^rnix). (4) 
1=1 i=l 

Hence, by applying the projector Pyv on both sides of (jl]) we further have 

M 

\fw) =^Ci\ui), where [wj) = Avl^i.m), and |/w) = Av|/v)- (5) 

i=l 

Denoting by Is the identity operator in S = V + W"*", the projector Pyy is obtainable from 
the relation Pyy = -^s — Pw^- Therefore the component |/>v) is available and can be used to 
determine the knots of the spline space to represent it. 

Let us suppose then that through the curvature function ([3]) we obtain a suitable partition 
for the space to represent l/w)- Now, in order to obtain the component |/v) from |/) in 
the most usual case involving subspaces V and W"*" close to each other, we need to find the 
representation of in a subspace of W = span{|Mj)}fl^. The approach for achieving such 
an aim is discussed in the next section. 



2.2 Determination of the signal representation through a nonUnear 
non-extensive approach 

At this point we can assume that we know the spanning set for W so that henceforth the 
problem is reduced to finding the representation with sparse coefficients. For this we could 
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apply the approach proposed in [4], which entails to use the normal equations 

M 

{Un\fw) = ^Ci{Un\Ui), n=l...,M. (6) 

1=1 

as constraints for the minimization of the norm like quantity ^fl^ IqI*^) < q < 1 but 
incorporating the equations in a stepwise manner. However, in the case motivating this work the 
number of necessary constraints is large enough to make the whole process slow. Hence rather 
than using the approach proposed in [4] we take an alternative route and apply a regularized 
version of the FOCUSS algorithm, which implies to minimize the functional 

M M 
>C = ^|c,r + A|||/n;)-^QK)|p, 

i=l i=l 

where A is a regularization parameter. The algorithm for implementing the approach is based 
on re- weighted least squares and is given in [22]. It comprises the following simple steps 

1) For each fixed q set a value for A and a value for the initial vector \c°) — Yl^i 

2) At each iteration, say iteration k, define the operators 

M M 
i=l i=l 

where the representation {x\Bi^m) i = 1 . . . , M of the kets i = 1 . . . , M are B-spline 

functions {x\Bi^rn) = Bi^rn{x) of order m. 

3) Compute Ic*^) as 

\^) = tiu{BlBu + \i)-^\fw), 
where 5^ indicates the adjoint of Bk and / the identity operator. 

4) Given a small e, while \\\c^) — \c^~^)\\ > e, repeat 2) and 3) 

For the derivation of the method and discussion on convergence issues see [22] . 

When the numerical convergence has been reached, say at iteration set q = cf , i — 
1, . . . , M and compute the required component |/v) as 

M 

l/v) ^^Ci\Bi^rn)- 
i=l 

3 Application to filtering of structured low frequency 
noise from a seismic signal 

We apply here the proposed approach to filtering low frequency noise from a seismic signal. As 
already mentioned, a common interference with broadband seismic signals is produced by long 

waves, generated by known or unknown sources, called infragravity waves [6]. This interference 
is referred to as low frequency noise as it falls in a frequency range of up to 0.05 Hz. Thus, the 



5 



model for the subspace of that type of structured noise, on a signal given by L = 403 samples, 
is 

= span{e^^=^^, ^ = 1, . . . , L}ll_,,. (7) 

The particular realization of the noise we have simulated is plotted in the top left graph of Fig 2 
(signal i = 1, . . . , L). However, it is appropriate to recall that the success of the approach 
does not depends on the actual form of the noise (as long as it belongs to the subspace given 
in ([7])) because the approach guarantees the suppression of the whole subspace W"*". 

The seismic signal ff,i = 1, . . . ,L shown in the right graph of Fig 2 is a piece of a test 
signal distributed by the seismic industry. The left bottom graph is the superposition of the 
signals in the top graphs, i.e. fi = fl^ + fi, i = I, . . . ,L. We fist subtract from / the component 
in W"*" to obtain /w = f — Pyv^f use this signal to find the points in the set ([3]). Then 
we subdivide uniformly those points to obtain 341 nonuniform knots defining the signal space. 
Using these knots we construct the nonuniform B-spline basis for the space (MATLAB codes 
for the implementation of both steps are available from [23]). 





Figure 2: Top left graph: Simulated low frequency noise (/"). Top right graph: piece of seismic 
signal distributed as test signal by the seismic industry (/''). Left bottom graph: Signal plus noise 
(/ = /* + /")■ Right bottom graph: Approximation recover from the left graph by applying the 
proposed approach for q = 0.123. 

The methodology discussed in the previous section requires to fix the values for q and A. To 
allow for good resolution the regularization parameter A is given a small value A = 10~^. As 
for the parameter q we have let it vary in the interval (0, 1] with step 0.001 and the best value 
resulted to be g = 0.123. The right bottom graph of Fig 2 depicts the filtered signal arising by 
applying the regularized FOCUSS method for A = 10"*^ and q = 0.123. 

The left graph of Fig 3 plots, |/'^ — the absolute value of the difference between the 
approximation (for q = 0.123) and the true signal For the sake of comparison, in the right 
graph we have plotted \ f^ — f^\, where is the approximation arising by filtering with Fast 
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Fourier Transform (FFT). For this approximation we simply take the FFT of /, ehminate the 
frequencies components in ([7]), and apply the inverse transform to obtain f^. The comparison 
shows the superiority of the proposed approach with respect to standard FFT filtering. 




0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



Figure 3: Left graph: Absolute value of the difference between the true seismic signal and the 
approximation obtained by the proposed approach with q = 0.123. Right graph: Absolute value 
of the difference between the true seismic signal and the approximation obtained by filtering 
frequencies in the FFT of /. 

In order to analyze the dependence of the solution on the parameter q we calculated the 
error's norm e(g) = ||/'^ — /*||. The plot of this error, against q, is depicted in Fig 4. While there 
is clearly an optimum value of q, the one we have used in this example (and a second best value 
q = 0.312) it should be mentioned that for all the values of q in (0 1] the approximation is 
superior than that obtained by filtering with FFT. 

1.4 1 , , , , , , , , , 1 
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Figure 4: Norm of the approximation error e{q) = \\fs — fq\\ as a function of the parameter q. 

As a final remark it may be worth stressing that, except for q = 1, the g-norm like quantity 
is not convex. Nevertheless, as opposed to the non-extensive g-entropy [9, 10] the g-norm like 
quantity is not extensive for g = 1. 

4 Conclusions 

The problem of structured noise suppression has been considered by modehing the subspaces 
of the signal components and applying a nonextensive nonlinear technique for separating them. 
The work was motivated by the problem of filtering infragravity waves from broadband seismic 
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signals. For this, the noise subspace was modelled using low frequency Fourier functions and 
that of the other component by a dedicated spline space (adapted to the signal in hand). A 
simulation involving a piece of seismic signal distributed by the seismic industry and noise up to 
0.05Hz has produced encouraging results in comparison with those arising by standard Fourier 
Transform filtering. 
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